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Abstract 

We analyse an epidemiological model of competing strains of pathogens and 
hence differences in transmission for first versus secondary infection due to inter- 
action of the strains with previously aquired immunities, as has been described 
for dengue fever (in dengue known as antibody dependent enhancement, ADE). 
Such models show a rich variety of dynamics through bifurcations up to deter- 
ministic chaos. Including temporary cross-immunity even enlarges the parameter 
range of such chaotic attractors, and also gives rise to various coexisting attrac- 
tors, which are difficult to identify by standard numerical bifurcation programs 
using continuation methods. A combination of techniques, including classical bi- 
furcation plots and Lyapunov exponent spectra has to be applied in comparison 
to get further insight into such dynamical structures. Here we present for the 
first time multi-parameter studies in a range of biologically plausible values for 
dengue. The multi-strain interaction with the immune system is expected to also 
have implications for the epidemiology of other diseases. 
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1 Introduction 

Epidemic models are classically phrased in ordinary differential equation (ODE) sys- 
tems for the host population divided in classes of susceptible individuals and infected 
ones (SIS system), or in addition, a class of recovered individuals due to immunity 
after an infection to the respective pathogen (SIR epidemics). The infection term in- 
cludes a product of two variables, hence a non-linearity which in extended systems 
can cause complicated dynamics. Though these simple SIS and SIR models only show 



fixed points as equilibrium solutions, they already show non-trivial equilibria arising 
from bifurcations, and in stochastic versions of the system critical fluctuations at the 
threshold. Further refinements of the SIR model in terms of external forcing or dis- 
tinction of infections with different strains of a pathogen, hence classes of infected with 
one or another strain recovered from one or another strain, infected with more than 
one strain etc., can induce more complicated dynamical attractors including equilibria, 
limit cycles, tori and chaotic attractors. 

Classical examples of chaos in epidemiological models are childhood diseases with 
extremely high infection rates, so that a moderate seasonal forcing can generate Feigen- 
baum sequences of period doubling bifurcations into chaos. The success in analysing 
childhood diseases in terms of modelling and data comparison lies in the fact that they 
are just childhood diseases with such high infectivity. Otherwise host populations can- 
not sustain the respective pathogens. In other infectious diseases much lower forces 
of infection have to be considered leading to further conceptual problems with noise 
affecting the system more than the deterministic part, leading even to critical fluctu- 
ations with power law behaviour, when considering evolutionary processes of harmless 
strains of pathogens versus occasional accidents of pathogenic mutants [1]. Only ex- 
plicitly stochastic models, of which the classical ODE models are mean field versions, 
can capture the fluctuations observed in time series data [2]. 

More recently it has been demonstrated that the interaction of various strains on 
the infection of the host with eventual cross-immunities or other interactions between 
host immune system and multiple strains can generate complicated dynamic attractors. 
A prime example is dengue fever. A first infection is often mild or even asymptomatic 
and leads to life long immunity against this strain. However, a subsequent infection 
with another strain of the virus often causes clinical complications up to life threaten- 
ing conditions and hospitalization, due to ADE. More on the biology of dengue and its 
consequences for the detailed epidemiological model structure can be found in Aguiar 
and Stollenwerk [3] including literature on previous modelling attempts, see also [4]. 
On the biological evidence for ADE see e.g. [5]. Besides the difference in the force of 
infection between primary and secondary infection, parametrized by a so called ADE 
parameter <j), which has been demonstrated to show chaotic attractors in a certain 
parameter region, another effect, the temporary cross-immunity after a first infection 
against all dengue virus strains, parametrized by the temporary cross- immunity rate 
q, shows bifurcations up to chaotic attractors in a much wider and biologically more 
realistic parameter region. The model presented in the Appendix has been described in 
detail in [3] and has recently been analysed for a parameter value of a = 2 year^ 1 cor- 
responding to on average half a year of temporary cross immunity which is biologically 
plausible [6] . For increasing ADE parameter <p first an equilibrium which bifurcates via 
a Hopf bifurcation into a stable limit cycle and then after further continuation the limit 
cycle becomes unstable in a torus bifurcation. This torus bifurcation can be located 
using numerical bifurcation software based on continuation methods tracking known 
equilibria or limit cycles up to bifurcation points [7]. The continuation techniques and 
the theory behind it are described e.g. in Kuznetsov [8]. Complementary methods 
like Lyapunov exponent spectra can also characterize chaotic attractor [9, 10], and led 



ultimately to the detection of coexisting attractors to the main limit cycles and tori 
originated from the analytically accessible fixed point for small <fi. Such coexisting 
structures are often missed in bifurcation analysis of higher dimensional dynamical sys- 
tems but are demonstrated to be crucial at times in understanding qualitatively the 
real world data, as for example demonstrated previously in a childhood disease study 
[11]. In such a study first the understanding of the deterministic system's attractor 
structure is needed, and then eventually the interplay between attractors mediated by 
population noise in the stochastic version of the system gives the full understanding of 
the data. Here we present for the first time extended results of the bifurcation struc- 
ture for various parameter values of the temporary cross immunity a in the region of 
biological relevance and multi-parameter bifurcation analysis. This reveals besides the 
torus bifurcation route to chaos also the classical Feigenbaum period doubling sequence 
and the origin of so called isola solutions. The symmetry of the different strains leads 
to symmerty breaking bifurcations of limit cycles, which are rarely described in the 
epidemiological literature but well known in the biochemical literature, e.g for cou- 
pled identical cells. The interplay between different numerical procedures and basic 
analytic insight in terms of symmetries help to understand the attractor structure of 
multi-strain interactions in the present case of dengue fever, and will contribute to the 
final understanding of dengue epidemiology including the observed fluctuations in real 
world data. In the literature the multi-strain interaction leading to deterministic chaos 
via ADE has been described previously, e.g. [12, 13] but neglecting temporary cross 
immunity and hence getting stuck in rather unbiological parameter regions, whereas 
more recently the first considerations of temporary cross immunity in rather compli- 
cated and up to now not in detail analysed models including all kinds of interations 
have appeared [14, 15], in this case failing to investigate closer the possible dynamical 
structures. 

2 Dynamical system 

The multistrain model under investigation can be given as an ODE system 



for the state vector of the epidemiological host classes x:= (S 1 , I\, I2, R) tr and besides 
other fixed parameters which are biologically undisputed the parameter vector of varied 
parameters a = (a,(p) tr . For a detailed description of the biological content of state 
variables and parameters see [3]. The ODE equations and fixed parameter values are 
given in the appendix. The equilibrium values x* are given by the equilibrium condition 
f{x*,a) = 0, respectively for limit cycles x*(t + T) = x*(t) with period T. For chaotic 
attractors the trajectory of the dynamical system reaches in the time limit of infinity 
the attractor trajectory x*(t), equally for tori with irrational winding ratios. In all 
cases the stability can be analysed considering small perturbations Ax(t) around the 
attractor trajectories 



— x = l(x,a) 



(1) 




• Ax 
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Here, any attractor is notified by x*(t), be it an equilibrium, periodic orbit or chaotic 
attractor. In this ODE system the linearized dynamics is given with the Jacobian 
matrix ^= of the ODE system Eq. (1) evaluated at the trajectory points x*(t) given 
in notation of [df /dx)\ xt . t y The Jacobian matrix is analyzed for equilibria in terms 
of eigenvalues to determine stability and the loss of it at bifurcation points, negative 
real part indicating stability. For the stability and loss of it for limit cylces Floquet 
multipliers are more common (essentially the exponentials of eigenvalues), multipliers 
inside the unit circle indicating stability, and where they leave eventually the unit circle 
determining the type of limit cycle bifurcations. And for chaotic systems Lyapunov 
exponents are determined from the Jacobian around the trajectory, positive largest 
exponents showing deterministic chaos, zero largest showing limit cycles including tori, 
largest smaller zero indicating fixed points. 



2.1 Symmetries 



To investigate the bifurcation structure of the system under investigation we first ob- 
serve the symmetries due to the multi-strain structure of the model. This becomes 
important for the time being for equilibria 1 and limit cycles. We introduce the follow- 
ing notation: With a symmetry transformation matrix S 
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we have the following symmetry: 
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with x* equilibrium values or x* = x*(t) limit cycle for all times t € [0, T\. For the 
right hand side / of the ODE system Eq. (1) the kind of symmetry found above is 



1 Equilibria are often called fixed points in dynamical systems theory, here we try to avoid this term, 
since in symmetry the term fixed is used in a more specific way, see below. 



called Z2-symmetry when the following equivariance condition holds 



f(Sx,a) = Sf(x,a) 



(5) 



with S a matrix that obeys S ^1 and S 2 = I, where I is the unit matrix. Observe that 
besides S also I satisfies (5). The symmetry transformation matrix S in Eq. (3) fulfills 
these requirements. It is easy to verify that the Z2-equivariance conditions Eq. (5) and 
the properties of S are satisfied for our ODE system. In Seydel [16] a simplified version 
of the famous Brusselator that shows this type of symmetry is discussed. There, an 
equilibrium and also a limit cycle show a pitchfork bifurcation with symmetry breaking. 

An equilibrium x* is called fixed when Sx* = x* (see [8]). Two equilibria x*,y* 
where Sx* ^ x*, are called S-conjugate if their corresponding solutions satisfy y* = Sx* 
(and because S 2 = I also x* = Sy*). For limit cycles a similar terminology is introduced. 
A periodic solution is called fixed when Sx*(i) = x*(t) and the associated limit cycles 
are also called fixed [8]. There is another type of periodic solution that is not fixed but 
called symmetric when 



where T is the period. Again the associated limit cycles are also called symmetric. 
Both types of limit cycles L are S-invariant as curves : SL = L. That is, in the 
phase-plane where time parameterizes the orbit, the cycle and the transformed cycle 
are equal. A S-invariant cycle is either fixed or symmetric. Two noninvariant limit 
cycles (SL ^ L) are called S-conjugate if their corresponding periodic solutions satisfy 
y*(t) = Sx*(t), Vi G M. The properties of the symmetric systems and the introduced 
terminology are used below with the interpretation of the numerical bifurcation analysis 
results. We refer to [8] for an overview of the possible bifurcations of equilibria and 
limit cycles of Z 2 -equivariant systems. 

3 Bifurcation diagrams for various a values 

We show the results of the bifurcation analysis in bifurcation diagrams for several 
a values, varying (ft continuously Besides the previously investigated case of a = 
2 year -1 , we show also a case of smaller and a case of larger a value, obtaining more 
information on the bifurcations possible in the model as a whole. The above mentioned 
symmetries help in understanding the present bifurcation structure. 

3.1 Bifurcation diagram for ct = 3 

For a = 3 the one-parameter bifurcation diagram is shown in Fig. 1 a). Starting 
with (ft = there is a stable fixed equilibrium, fixed in the above mentioned notion 
for symmetric systems. This equilibrium becomes unstable at a Hopf bifurcation H at 
(ft = 0.164454. A stable fixed limit cycle originates at this Hopf bifurcation. This limit 
cycle shows a supercritical pitch- fork bifurcation P~ , i.e. a bifurcation of a limit cycle 
with Floquet multiplier 1, splitting the original limit cycle into two new ones. Besides 




(6) 



the now unstable branch two new branches originate for the pair of conjugated limit 
cycles. The branches merge again at another supercritical pitch-fork bifurcation P~, 
after which the limit cycle is stable again for higher 0- values. The pair of S-conjugate 
limit cycles become unstable at a torus bifurcation TR at (p = 0.89539. 
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Figure 1: a) a — 3: Equilibria or extremum values for limit cycles for logarithm of total infected 
h + h + I\2 + hi- Solid lines denote stable equilibria or limit cycles, dashed lines unstable 
equilibria or periodic-one limit cycles. Hopf bifurcation H around <p — 0.16 two pitchfork 
bifurcations P~ and a torus bifurcation TR. Besides this main bifurcation structure we found 
coexisting tangent bifurcations T between which some of the isolas live, see especially the one 
between <f) = 0.71 and 0.79. Additionally found flip bifurcations are not marked here, see text, 
b) a = 2: In this case we have a Hopf bifurcation H at <f> = 0.11, and besides the similar 
structure as found in a) also more separated tangent bifurcations T at (j> = 0.494, 0.539, 0.931, 
0.978 and 1.052 c) a = 1: Here we have the Hopf bifurcation at <fi = 0.0598 and thereafter 
many tangent bifurcations T, again with coexisting limit cylces. 



Besides this main bifurcation pattern we found two isolas, that is an isolated solu- 
tion branch of limit cycles [17]. These isola cycles L are not S-invariant, that is SL ^ L. 
Isolas consisting of isolated limit cycles exist between two tangent bifurcations. One 
isola consists of a stable and an unstable branch. The other shows more complex bi- 
furcation patterns. There is no full stable branch. For (f> = 0.60809 at the tangent 
bifurcation T a stable and an unstable limit cycle collide. The stable branch becomes 
unstable via a flip bifurcation or periodic doubling bifurcation F, with Floquet multi- 
plier (—1), at 4> = 0.61918 which is also pitchfork bifurcation for the period-two limit 
cycles. At the other end of that branch at the tangent bifurcation T at <f> = 0.89768 
both colliding limit cycles are unstable. Close to this point at one branch there is a 
torus bifurcation TR, also called Neimark-S acker bifurcation, at <f> = 0.89539 and a flip 
bifurcation F at cf) = 0.87897 which is again a pitchfork bifurcation P for the period- 
two limit cycles. Contiuation of the stable branch originating for the flip bifurcation F 
at (f) = 0.61918 gives another flip bifurcation F at cf> = 0.62070 and one closed to the 
other end at cp = 0.87897, namely at (f> = 0.87734. These results suggest that for this 
isola two classical routes to chaos can exist, namely via the torus or Neimark-Sacker 
bifurcation where the dynamics on the originating torus is chaotic, and the cascade of 
period doubling route to chaos. 



3.2 Bifurcation diagram for a = 2 



For a = 2 the one-parameter bifurcation diagram is shown in Fig. 1 b) . The stable fixed 
equilibrium becomes unstable at a supercritical Hopf bifurcation H at <p = 0.1132861 
where a stable fixed limit cycle originates. This stable limit cycle becomes unstable 
at a superciritcal pitchfork bifurcation point P~ at (f> = 0.4114478 for a limit cycle. 
This point marks the origin of a pair of S-conjugate stable limit cycles besides the now 
unstable fixed limit cycle. Here one has to consider the two infected subpopulations 
I\ and I2 to distinguish the conjugate limit cycles. Because the two variables I\ and 
I2 are interchangeable this can also be interpreted as the stable limit cycles for the 
single variable say I\ . The fixed stable equilibrium below the Hopf bifurcation where 
we have I* = 7|, R\ = R 2 , S\ = S 2 and I* 2 = I 21 is a fixed equilibrium. For 
the fixed limit cycle in the parameter interval between the Hopf bifurcation and the 
pitchfork bifurcation we have Jf(t) = J|(t), R{(t) = R^t), S$(t) = S%(t) and I{ 2 {t) = 
/|i(t). This means that at the Hopf bifurcation H the stable fixed equilibrium becomes 
an unstable fixed equilibrium. In the parameter interval between the two pitchfork 
bifurcations P~ at (f> = 0.4114478 and subcritical P + at <j) = 0.9921416, two stable limit 
cycles coexist and these limit cycles are S-conjugate. At the pitchfork bifurcation points 
the fixed limit cycle becomes unstable and remains fixed, and two stable S-conjugate 
limit cycles originate (see [8, Theorem 7.7]). The invariant plane I\ = I2, Ri = R2, S\ = 
S21I12 = I21 forms the separatrix between the pair of stable S-conjugate limit cycles 
x*(t) and Sx*(t), Vt € BL The initial values of the two state variables S(to) and R(to) 
together with the point on the invariant plane, determine to which limit cycle the system 
converges. Continuation of the stable symmetric limit cycle gives a torus or Neimark- 
Sacker bifurcation at point denoted by TR at <p = 0.5506880. At his point the limit 
cycles become unstable because a pair of complex-conjugate multipliers crosses the unit 
circle. Observe that at this point in the time series plot [3, there Fig. 12] the chaotic 
region starts. In [18] the following route to chaos, namely the sequence of Neimark- 
Sacker bifurcations into chaos, is mentioned. Increasing the bifurcation parameter <p 
along the now unstable pair of S-conjugate limit cycles leads to a tangent bifurcation 
T at (ft = 1.052418 where a pair of two unstable limit cycles collide. This branch 
terminates at the second pitchfork bifurcation point denoted by P + at <f> = 0.9921416. 
Because the first fold point gave rise to a stable limit cycle and this fold point to an 
unstable limit cycle we call the first pitchfork bifurcation supercritical and the latter 
pitchfork bifurcation subcritical. These results agree very well with the simulation 
results shown in the bifurcation diagram for the maxima and minima of the overall 
infected [3, there Fig. 15]. Notice that AUTO [7] calculates only the global extrema 
during a cycle, not the local extrema. Fig. 1 b) shows also two isolas similar to those 
for a = 3 in Fig. 1 a). 

3.3 Bifurcation diagram for a = 1 

For a = 1 the bifurcation diagram is shown in Fig 1 c) . In the lower <ft parameter range 
there is bistability of two limit cycles in an interval bounded by two tangent bifurcations 



$F,TR,T.P$ 




\ 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 i \ 3.75 3.77 3.70 3.81 3.83 3.85 \ 3.75 3.77 3.70 3.81 3.83 3.85 

a) ♦ bj s cj s 

Figure 2: a) a. = 1. Detail of Fig. 1 c). We find pitchfork bifurcations P at <f) — 0.239 and 
0.325, flip bifurcations F at <f> = 0.298, 0.328,0.344,0.346, 0.406, 0.407, 0.411 and 0.422, further 
tangent bifurcations T at <f> = 0.292, 0.346 and 0.422. Four almost coexisting bifurcations, 
namely F's at <j> = 0.4112590. b) and c) state space-plots of susceptibles and logarithm of 
infected for a = 1 and 4> = 0.294 for two coexisting stable limit cycles. 

T. The stable manifold of the intermediate saddle limit cycle acts as a separatrix. 
Inceasing <f> the stable limit cycles become unstable at the pitchfork bifurcation P at 
(j) = 0.2390695. Following the unstable primary branch, for larger values of 4> we observe 
an open loop bounded by two tangent bifurcations T. The extreme value for <p is at 
4> = 0.6279042. Then lowering (j) there is a pitchfork bifurcation P at <fi = 0.5016112. 
Later we will return to the description of this point. Lowering further the limit cycle 
becomes stable again at the tangent bifurcations T at 4> = 0.3086299. Increasing (j) this 
limit cycle becomes unstable again at the pitchfork bifurcation P at 4> = 0.3253242. 

Continuation of the secondary branch of the two S-conjugated limit cycles from 
this point reveals that the stable limit cycle becomes unstable at a torus bifurcation 
TR at cf) = 0.4257346. The simulation results depicted in [3, Fig. 13] show that there is 
chaos beyond this point. The secondary pair of S-conjugate limit cycles that originate 
from pitchfork bifurcation P at <fi = 0.2390695 becomes unstable at a flip bifurcation 
F. Increasing (ft further it becomes stable again at a flip bifurcation F. Below we 
return to the interval between these two flip bifurcations. The stable part becomes 
unstable at a tangent bifurcation T, then continuing, after a tangent bifurcation T 
and a Neimark-S acker bifurcation TR. This bifurcation can lead to a sequence of 
Neimark-S acker bifurcations into chaos. The unstable limit cycles terminates via a 
tangent bifurcation F where the primary limit cycle possesses a pitchfork bifurcation 
P at <p = 0.5016112. At the flip bifurcation F the cycle becomes unstable and a new 
stable limit cycle with double period emanates. The stable branch becomes unstable at 
a flip bifurcation again. We conclude that there is a cascade of period doubling route 
to chaos. Similarly this happens in reversed order ending at the flip bifurcation where 
the secondary branch becomes stable again. 

Fig. 2 a) gives the results for the interval 0.28 < 4> < 0.44 where only the minima 
are show. In this plot also a "period three" limit cycle is shown. In a small region it 
is stable and coexists together with the "period one" limit cycle. The cycles are shown 
in Fig. 2 b) and c) for <fi = 0.294. The one in c) looks like a period-3 limit cycle. In 
Fig. 2 continuation of the limit cycle gives a closed graph bounded at the two ends by 
trangent bifurcations T where a stable and an unstable limit cycle collide. The intervals 



where the limit cycle is stable, are on the other end bounded by flip bifurcations F. 
One unstable part intersects the higher period cycles that originate via the cascade of 
period doubling between the period-1 limit cycle flip bifurcations F at (ft = 0.3281636 
and (ft = 0.4112590. This suggest that the period-3 limit cycle is associated with a 
"period-3 window" of the chaotic attractor. We conjecture that this interval is bounded 
by two homoclinic bifurcations for a period-3 limit cycle (see [19, 20, 21, 22]). The 
bifurcation diagram shown in [3, there Fig. 13] shows the point where the chaotic 
attractor disappears abruptly, possible at one of the two homoclinic bifurcations. In 
that region the two conjugated limit cycles that originate at the pitchfork bifurcation 
P at (ft = 0.3253242 are the attractors. These results suggest that there are chaotic 
attractors associated with the period-1 limit cycle, one occurs via a cascade of flip 
bifurcations originating from the two ends at (ft = 0.3281636 and (ft = 0.4112590 and 
one via a Neimark-Sacker bifurcation TR at (ft = 0.4257346. 



4 Two- parameter diagram 

We will now link the three studies of the different a values by investigating a two- 
parameter diagram for (ft and a, concentrating especially on the creation of isolated 
limit cycles, which sometimes lead to further bifurcations inside the isola region. Fig. 3 
gives a two-parameter bifurcation diagram where (ft and a are the free parameters. For 
low (ft- values there is the Hopf bifurcation H and all other curves are tangent bifurcation 
curves. 
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Figure 3: Two-dimensional parameter bifurcation diagram with (ft and a as parameters. Only 
one Hopf bifurcation (dotted lines) and many tangent bifurcation curves (dashed lines) are 
shown in the range a G [1,4]. The isolated limit cycles originate above a = 3. For lower values 
of a periodic doubling routes to chaos originate. 



Isolas appear or disappears upon crossing an isola variety. At an elliptic isola 
point an isolated solution branch is born, while at a hyperbolic isola point an isolated 
solution branch vanishes by coalescence with another branch [17]. From Fig. 3 we see 
that at two values of a > 3 isolas are born. Furthermore, period doubling bifurcations 
appear for lower a values, indicating the Feigenbaum route to chaos. However, only 
the calculation of Lyapunov exponents, which are discussed in the next section, can 
clearly indicate chaos. 



5 Lyapunov spectra for various at values 



The Lyapunov exponents are the logarithms of the eigenvalues of the Jacobian matrix 
along the integrated trajectories, Eq. (2), in the limit of large integration times. Besides 
for very simple iterated maps no analytic expressions for chaotic systems can be given 
for the Lyapunov exponents. For the calculation of the iterated Jacobian matrix and 
its eigenvalues, we use the QR decomposition algorithm [23]. 




a) ♦ b) * c) * 

Figure 4: Spectrum of the four largest Lyapunov exponents with changing parameter <j) and a) 
fixed a = A, b) a = 2 and c) a = 1. 

In Fig. 4 we show for various a values the four largest Lyapunov exponents in 
the (f) range between zero and one. For a = 4 in Fig. 4 a) we see for small <f> values 
fixed point behaviour indicated by a negative largest Lyapunov exponent up to around 
4> = 0.2. There, at the Hopf bifurcation point, the largest Lyapunov exponent becomes 
zero, indicating limit cycle behaviour for the whole range of (f>, apart from the final bit 
before <f> = 1, where a small spike with positive Lyapunov exponent might be present, 
but difficult to distinguish from the noisy numerical background. 

For a = 2 in Fig. 4 b) however, we see a large window with positive largest 
Lyapunov exponent, well separated from the second largest being zero. This is s clear 
sign of deterministically chaotic attractors present for this <f> range. Just a few windows 
with periodic attractors, indicated by the zero largest Lyapunov exponent are visible 
in the region of 0.5 < <j> < 1. For smaller <fi values we observe qualitatively the same 
behaviour as already seen for a = 4. For the smaller value of a = 1 in Fig. 4 c) the 
chaotic window is even larger than for a = 2. Hence deterministic chaos is present for 
temporary cross immunity in the range around a = 2 year -1 in the range of <f> between 
zero and one. 



6 Conclusions 

We have presented a detailed bifurcation analysis for a multi-strain dengue fever model 
in terms of the ADE parameter <j), in the previously not well investigated region between 
zero and one, and a parameter for the temporary cross immunity a. The symmetries 
implied by the strain structure, are taken into account in the analysis. Many of the 
possible bifurcations of equilibria and limit cycles of Z2-equivariant systems can be 
distinguished. Using AUTO [7] the different dynamical structures were calculated. 



Future time series analysis of epidemiological data has good chances to give insight 
into the relevant parameter values purely on topological information of the dynamics, 
rather than classical parameter estimation of which application is in general restricted 
to farely simple dynamical scenarios. 
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7 Appendix: Epidemic model equations 



The complete system of ordinary differential equations for a two strain epidemiological 
system allowing for differences in primary versus secondary infection and temporary 
cross immunity is given by 



For two different strains, 1 and 2, we label the SIR classes for the hosts that have 
seen the individual strains. Susceptibles to both strains (S) get infected with strain 1 
(ii) or strain 2 (I2), with infection rate (3. They recover from infection with strain 1 
(becoming temporary cross- immune R\) or from strain 2 (becoming R2), with recovery 
rate 7 etc.. With rate a, the Ri and R 2 enter again in the susceptible classes (Si 
being immune against strain 1 but susceptible to 2, respectively 52), where the index 
represents the first infection strain. Now, Si can be reinfected with strain 2 (becoming 
/12), meeting I2 with infection rate (5 or meeting I12 with infection rate (p(3, secondary 
infected contributing differently to the force of infection than primary infected, etc.. 

We include demography of the host population denoting the birth and death rate 
by (i. For constant population size N we have for the immune to all strains R = 
N — (S + Ii + 12 + R\ + R2 + Si + S2 + 112 + hi) and therefore we only need to consider 
the first 9 equations of Eq. (7), giving 9 Lyapunov exponents. In our numerical studies 
we take the population size equal to N = 100 so that numbers of susceptibles, infected 
etc. are given in percentage. As fixed parameter values we take fi = (1/65) year -1 , 
7 = 52 year" 1 , j3 = 2 ■ 7. The parameters <\> and a are varied. 





7/2 - (a + n)R2 



(7) 




